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Abstract 

We develop 3rd order maximum-principle-satisfying direct discontinuous Galerkin methods 
0 eu na mi for convection diffusion equations on unstructured triangular mesh. We carefully 
calculate the normal derivative numerical flux across element edges and prove that, with proper 
choice of parameter pair (/?□, /3i) in the numerical flux, the quadratic polynomial solution satisfies 
strict maximum principle. The polynomial solution is bounded within the given range and third 
order accuracy is maintained. There is no geometric restriction on the meshes and obtuse 
triangles are allowed in the partition. A sequence of numerical examples are carried out to 
demonstrate the accuracy and capability of the maximum-principle-satisfying limiter. 

Keywords: Discontinuous Galerkin methods; Convection diffusion equation; Maximum Prin¬ 
ciple; Positivity Preserving; Incompressible Navier-Stokes equations; 


1 Introduction 

In this article, we study direct discontinuous Galerkin finite element method [S] and its variations 
mmm to solve two-dimensional convection diffusion equations of the form, 

ut + V • F(u) — V • (A(u)Vu) = 0, (x, y, t) e D x (0, T), (1.1) 

with zero or periodic boundary conditions. We have spacial domain HcK 2 and initial condition 
u(x,y, 0) = uo(x,y). The convection flux is denoted as F(u) = ( f(u),g(u )) and diffusion matrix 
A(u) = ( a,ij(u )) is assumed symmetric and positive definite. 

On the continuous level, solution of 0 may satisfy the maximum principle, which states 
the evolution solution u(x, y, t ) being bounded below and above by the given constants, m < 
u(x, y, t) < M. Here m and M are the lower and upper bounds of the initial and boundary data. 
It is desirable that the numerical solution satisfies the discrete maximum principle. The discrete 
maximum principle can be considered as a strong L°° sense stability result. Failure of preserving 
the bounds or maintaining the positivity of the numerical solution may lead to ill-posed problems 
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and practically cause the computations to blow up. Thus it is attractive to have the numerical 
solution satisfy discrete maximum principle (or preserve positivity). Solution of equation 0 
may represent a specific physical meaning and is supposed to be positive, thus negative value 
approximation loses physical meanings in such cases. 

Generally it is very difficult to design high order numerical methods that satisfy discrete max¬ 
imum principle for convection diffusion equations (1.1). No finite difference method is known to 
achieve better than second-order accuracy 0ES1I2O] that satisfies discrete maximum principle. 
Much less is known for higher order methods such as spectral FEM, hp-FEM or finite volume 
methods 0 e na Eg. Compared to the elliptic type, more restrictive conditions on mesh are 
required to obtain discrete maximum principle for the parabolic type equations, see usmsuzun]. 

In this article, we study direct discontinuous Galerkin method [5] and its variations mmm, 
and prove the polynomial solution satisfy discrete maximum principle with third order of accu¬ 
racy. Discontinuous Galerkin (DG) method is a class of finite element method that use completely 
discontinuous piecewise functions as numerical approximations. Since the basis functions can be 
completely discontinuous, these methods have the flexibility that is not shared by standard finite 
element methods, such as the allowance of arbitrary triangulations with hanging nodes, complete 
freedom of choosing polynomial degrees in each element (p adaptivity), and extremely local data 
structure and the resulting high parallel efficiency. 

Recently in [241 [25l [26], Zhang and Shu designed a maximum-principle-satisfying limiter for 
high order DG and finite volume methods for hyperbolic conservation laws. The key step in Zhang 
and Shu’s discussion is to show the polynomial solution average falling in the given minimum and 
maximum bounds. For hyperbolic type equations, the solution average evolution only relies on 
the solution polynomial values on the element edges. For diffusion type equations, the evolution 
of solution average depends on the solution derivative values on the edges, thus the technique 
developed in [23] can not be applied. 

In [8], we developed the direct DG method (DDG) as a new diffusion solver. The key con¬ 
tribution of direct DG method is the introduction of numerical flux to approximate the solution 
derivative at the discontinuous element boundaries. The scheme is directly based on the weak 
formulation of diffusion equation, thus gains its name the direct DG method. Now let’s use the 
simple 2-D heat equation to go through the main idea of direct DG method, 


ut — A u = 0. (1.2) 

Multiply the heat equation with test function v, integrate over the element K, have integration by 
parts and formally we obtain, 


utv dxdy — / u^v ds+ Vm • Vu dxdy = 0. 


IK 


' dK 


IK 


The numerical flux fqj introduced in [8] is defined as follows, 


_—~— I'll] 

tin = Vu • n = /5 0 y— 1 + 
h K 

It involves the jump, the derivative average and higher order derivative jumps to approximate the 
normal derivative u n on the element boundary dK. Here n = ( 711 , 712 ) is the outward unit normal 
along 3K and h k is the diameter of element K. The coefficient pair (/3q, /3i) is chosen to guarantee 
the convergence of the scheme. 


du 

T Pl O K [TtnnJ • 
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Due to accuracy loss of the original DDG method [8], we further developed DDG method with 
interface correction in [9] in which optimal (k + l)th order convergence is obtained with any order 
P k polynomial approximations. We also have the symmetric [19] and nonsymmetric version [21] 
of the DDG methods. In this paper, we mainly carry out the maximum principle study on DDG 
method with interface correction [9] since it is the most efficient solver for time dependent diffusion 
equations. The maximum principle arguments discussed in the following sections also apply to 
DDG method [E] and its symmetric and nonsymmetric variations H3EU. 

In [22], we prove the DDG solutions satisfy discrete maximum principle on rectangular and 
uniform triangular meshes with 3rd order of accuracy. We use an algebraic methodology and a 
monotonicity argument to show the polynomial solution average being bounded within the given 
range. The DG polynomial solution was written out in the Lagrange format with the unknowns 
carefully chosen on the element. With Euler forward in time, we show the solution average at next 
time level depends on the current time level solution values in a monotone fashion. For unstructured 
mesh with possible obtuse triangles, it is very hard to identify six degrees of freedom to represent 
the P 2 quadratic polynomial solution such that the monotone argument in [22j can be applied. 


In this article we extend maximum principle studies of (1.1) on unstructured triangle mesh 


Again let’s use the heat equation to illustrate the new technique to carry out the proof. Notice that 
the key step of the discussion is to show the solution average falling in the given range. Take test 
function v = 1 in the DDG scheme and discretize in time with Euler forward, we have the solution 
average evolving in time as, 


wjTi +1 — n 


U 


K 


= U'k + 


At 


area(K) J dK 


ds, 


with the average defined as v n K = are ^ K j f K u\(x, y ) dxdy and n^(x, y) as the polynomial solution 
at time step t n in element K. 

Instead of identifying suitable locations as degrees of freedom and writing out u^(x,y) in the 
Lagrange format as in [22], we directly calculate the normal derivative flux from the given 
solution values in element K and its neighbors. Given suitable choice of coefficient pair (/3o,/3i) 
in the numerical flux, we can bound the solution average E [m, M] at time level f n +i once 

we know u^(x,y) E [rn, M] at previous time level t n . Finally we borrow the maximum principle 


discussion of [25] to show the DG polynomial solution of general convection diffusion equations (1.1) 
satisfy strict maximum principle with 3rd order of accuracy. A sequence of numerical examples are 
carried out to demonstrate the DG solutions are strictly bounded by the given values and at the 
same time maintain the 3rd order accuracy. Solutions to nonlinear porous medium equations with 
nonnegative initial data are maintained sharply nonnegative. Examples of incompressible Navier- 
Stokes equations with high Reynold numbers are tested. Overshoot and undershoots are removed 
with maximum principle limiter applied. 

The key feature of direct DG methods is the introduction of numerical flux that approximates 
the solution derivative u n on the discontinuous element boundary. This gives direct DG methods 
the extra flexibility and advantage over IPDG method [T] and LDG method [1J. Following this 
maximum principle framework, both IPDG method and LDG method can be proved to satisfy 
maximum principle with up to second order of accuracy. In E3, DG solutions with piecewise linear 
polynomial approximations are shown satisfying maximum principle on unstructured triangle mesh. 

The paper is organized as follows. We first review the scheme formulation of direct DG method 
with interface correction [9] in section [2j In Section [3j we prove the direct DG solutions satisfy 
discrete maximum principle with 3rd order accuracy. We conduct numerical tests to validate the 
theoretical results in Section [4j Section [5] serves as the Appendix in which we provide one way to 
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construct a specific quadrature rule for quadratic polynomials with selected points as quadrature 
points. 


2 Direct DG method with interface correction 


We first recall the scheme formulation of direct DG method with interface correction [9] for two- 
dimensional diffusion equations, 


itt — V • (A(u)X7u) = 0, (x, y,t) G Cl x (0, T), 


( 2 . 1 ) 


with initial condition u(x, y, 0) = uq(x, y ) and zero or periodic boundary conditions. The complete 


scheme formulation of convection diffusion equation (1.1) will be laid out toward the end of this 


section. We should specify the DG method is for spatial discretization and we will incorporate 
high order TVD Runge-Kutta methods na nn to march forward the solution in time. As an 
explicit scheme, our method is thus more efficient for convection dominated problems. However, 
the extremely local dependency allows a very efficient parallelization and dramatically improves 
the efficiency of the explicit method. 

Let 7ft be a shape-regular partition of the polygonal domain Cl into triangle elements {K\ xeT fl 
with Cl = UK&T h K- By hx = diam(/\), we denote the diameter of the triangle element K G 7/j. 
We denote h = max^g^ hx as the mesh size of the partition. We have P k (K ) representing the 
A'th degree polynomial space on element K. The DG solution space is defined as, 

Y k h = {ve L\Cl ) : v\ K G P k (K),\/K G T h }. 

Suppose K and K' are two adjacent triangles and share one common edge e. There are two traces 
of v along the edge e, where we add or subtract those values to obtain the average and the jump. 
The outward normal vector from K to its neighbor element K' is denoted by n = (711,712). Now 
the average and jump of v on the edge e are defined as follows, 

v = ^ (v\ K + v\x '), [v\ = v\ K > - v\ K . 


The original DDG scheme of (2.1) defined in [8] is to find DG solution u G Y k , such that for 


any test function v G Yf we have, 


/ utv dxdy — / (A(u)Vrt • n)u ds + / A(u)Yu ■ Yv dxdy = 0, V K G Th- (2.2) 
Ik JdK Jk 


The numerical flux A(u)Yu ■ n (equation (3.7) of |21] in dimension-by-dimension format) along the 
element edge is defined as, 


A(u)Yu ■ n = \ bn(u) x + b l 2 (u) y j m + [b 2 i(u) x + b 2 2{u) y j n 2 , 

where bij(u ) = / a l3 (u)du with aij{u ) as the diffusion matrix A{u) entry. The outward normal is 
given with n = (ni, 772 ). Similar to (3.7) of [21], for example, the numerical flux b\\(u) x is calculated 
with formula, 

bu(u) x = + bn(u) x + fah K {[b\i{u) xx ]ni + [bn(u) yx \n 2 } . 

n K 
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The coefficient pair (/3o,/?i) should be chosen carefully to ensure the stability and convergence of 
the scheme. Again we have Hr as element IPs diameter or the length of edge dK. Notice in 
we further require the entries aij{u) > 0 to carry out the discrete maximum principle proof. 

In this article, we simplify the calculation of numerical flux A(u)X7u ■ n to the following, 


u 


A(u)Vu • n = Vu • 7 = (3q- - 1 - rt 7 + fdihxlu^A, 

Hr 


(2.3) 


where 7 = A T (u )n is a vector pointing from element K into its neighbor along the edge dK. The 
simplification holds true since (A(u)Vtt) • n = Vrt • (A T (rt)n) and A{u) is positive definite. We 
should point out the simplified version (|2.3|) is very important to carry out the maximum principle 


discussion in the following sections. Now the DDG interface correction [9] of (2.1) is defined to find 
solution u 6 V^, such that for any test function t; £ ¥j we have, 


utv dxdy — 


IK 


/ (A('u)Vw ■ n)w ds + / A(u)Vu ■ V v dxdy + / A(v)Vv ■ n[it] ds = 0, V K e Th, 

IdK J K JdK 

(2.4) 


with the numerical flux (2.3). The last term in (2.4) is the the extra added interface correction. 


Notice that the test function v is taken to be zero outside the element K, thus the derivative average 
degenerates to A(v)Vv = ^ A(v)Vv\r on the edge dK. 


The complete scheme formulation of (1.1) with DDG interface correction follows, 


/ UfV dxdy+ / F ■ nv ds — F ■ Vu dxdy 

Ik JdK Jk 


= / (A(u)Vu ■ n)v ds — A(u)Vu ■ Vu dxdy — / A(u)Vu • n[u] ds, (2.5) 
JdK Jk JdK 

with the convection term Lax-Friedrichs flux defined as, 

F ■ n = — (F(ur) ■ n + F(ur/) ■ n — ol(uri — ur)) , with a = max |F'(u) • n|. 

2 U 

3 Maximum-Principle-Satisfying DDG methods 


In this section we prove DDG polynomial solutions of nonlinear diffusion equations (2.1) satisfy 
discrete maximum principle with the M-P-S limiter applied. We first discuss the linear case on 


unstructured triangle mesh in section 3.1 Then we extend the study to nonlinear diffusion equations 
in section 13.21 

Notice that the second derivative jump term has no contribution to the calculation of the 
numerical flux with low order P° and P 1 approximations. The scheme of DDG method with 
interface correction (2.4) degenerates to IPDG methods with low order approximations. In this 


paper, we focus on P^quadratic polynomial approximations with 3rd order of accuracy. We skip 
the trivial piecewise constant case and refer to m for 2 rd order linear approximations. 

On the continuous level the maximum principle states that m < u{x, y, t ) < M , given m and 
M as the minimum and maximum of the initial data u(x,y, 0) = uo{x,y ) and boundary data. 
We have u(x,y,t n ) to denote the exact solution at time level t n and u^(x, y) as the polynomial 
solution on element I\ and at time level t n . Our goal is to prove the polynomial solution satisfy 
m < Uf((x,y) < M without losing the 3rd order accuracy at all time levels. We can simplify the 
discussion to Euler forward time discretization, since the full scheme (high order strong stability 
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preserving (SSP) Runge-Kutta method) is a convex combination of Euler forward scheme. For 
example, the third order SSP Runge-Kutta method in m is 

f =u n + AtH(u n ) 

< = \u n + + A tH(u^)) 

{ u n+1 = \u n + l( U W + A tH( U W)) 

Now assume at time level t n , we have 1) the DDG solution is 3rd order accurate; and 2) we 
have u v K {x. y) G [rri. M] on all elements. The goal is to prove the solution polynomial u r £ rl (x, y ) at 
next time level t n +i still stay inside the bounds [m, M] without losing accuracy. To carry out this 
study we need to consider following two steps: 

1 . to prove the polynomial solution average u 1 ^ 1 stay inside the bounds [m,M ]; 

2. to prove the whole polynomial u r ^' 1 (x,y) stay inside [m,M] without losing accuracy. 

The most challenging and the major step is to show the polynomial average u falling in 
[m, M]. For the second step, we simply apply a linear scaling limiter [11] to y ) and obtain a 

modified polynomial n , ^ 1 (x,y) such that the whole polynomial u’^ 1 (x,y) G [ m,M] without losing 
accuracy, we refer to [25] for the proof of this accuracy preserving limiter. We have the DDG 
solution approximates the exact solution with 3rd order accuracy, thus the polynomial solution can 
only jump out of the bounds [m, M] in the scale of h 3 with h as the mesh size. The limiter m is 
applied to compress and squeeze the polynomial in the scale of h 3 and put it back into the bounds 
[m, M]. The extra cost to preserve maximum principle is to apply the limiter illj to maintain the 
bounds. 


3.1 Linear diffusion equation 


In this section, we prove the DDG quadratic polynomial solutions of heat equation (1.2) satisfy 
discrete maximum principle. Again we focus on the first step and investigate under what conditions 
the solution average u 1 ^ 1 G [m,M] given u^(x,y) G [m,M] on all e lements. 

The scheme formulation of DDG with interface correction for (1.2) is to find DG solution u G 
such that for any test function v G we have, 


/ utv dxdy — / u^v ds + / V u ■ V v dxdy + / u n [tt] ds = 0, VA'G7h. (3.1) 
IK JdK JK JdK 


The numerical flux Iqj on the element boundary dK is given with, 


^^ yli\ 

U n = Pq— b u n + P\ hx [finn] • 

h K 

To obtain the solution average evolution, we take test function v 
with forward Euler and formally we have, 


—n+l _ — n 


U 


K 


~ U K + 


At 


area{K ) J dK 


ds. 


(3.2) 


1 in (3.1), discretize in time 


(3.3) 


The solution average is given with u\ = are ^ K ^ Jx u k( x i d) dxdy and At is denoted as the time 
step size. The average can be calculated out exactly, since u'j ( (x,y) is a quadratic polynomial. 
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Thus we see the quantity uff 
edges. Recall that the numerical flux u 


:r ' n+1 of (3.3) is essentially determined by the integral of u n on the three 
of (3.2) involves the solution jump, normal derivative 

eventually is a 


— n+1 


average and second order normal derivative jumps on dK. The quantity u K 
function of the four solution polynomials that spread out in K and its three neighbors. 

For uniform triangular mesh in [22], we pick six solution values on each element and write 
out the P 2 solution polynomial in Lagrange format, then use a monotone argument to bound 
E [m,M ]. For arbitrary triangular mesh, it’s hard to identify such six points inside each 
triangle. We will use a new idea to bound 1 . 

We observe that the three quantities, namely u, u n and u nn restricted on the edges from each 
element, contribute to the calculation of u^. If we manage to calculate u, u n and u nn from the 
given solution polynomials, we can easily bound u ^ +1 E [m,M] once we have u 1 f(x,y) E [m,M] 
for all K. 


C 




Figure 1: Left: K and its neighbor elements. Right: selected points to calculate rqj on edge AB. 


Let’s use Figure [l] to illustrate the solution points selected to calculate u. u n and u nn . For 
example, we consider the edge AB shared by I\ and its neighbor element K 3 (the right one in 
Figure [lj). Notice that u r f 3 (x,y) is a P 2 polynomial. The three points, namely x^ 3 i i, x^ 3 j 2 , x^ 3 i 3 , 
are enough to represent the restriction of u r f 3 (x. y) on edge AB. The second normal derivative 


Unn degenerates to a constant, and we can use three points x'^ 3 , x^- 34 , 


x ^ 3 5 on the normal 


line through the edge center to calculate u n n- The first normal derivative u n on AB is a linear 
polynomial, thus the same three points on the normal line are good enough to calculate the line 
integral of u n - With this new idea to calculate the numerical flux ujj of (3.3), we are ready to 
bound the average u ? +1 


A K 


Theorem 3 . 1 . Consider DDG scheme with interface correction (3.1) - (3.2) with P 2 quadratic 
approximations on unstructured triangular mesh. Given u 7 f i (x,y) in the range of[m,M] for all K, 


we have uf ^ 1 E [ m,M] provided, 


Po>l~6(h, g < A < \- 


A = 


At 


aread(K) 


<A{p 0 ,A,M). 


(3.4) 


Here (/3q,/5i) is the coefficient pair in the numerical flux (3.2). We have 6 and 0 denoted as the 
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maximum and minimum angle of the partition Th, and A is a function of /do, f}\, 0 and 0, i.e. 


A = tan(0) • min ■ 


W\ 


tan( 0 ) 


mm 


w 1 


w\ 


72(1 - 4/30 ’ tan( 0 ) l 6(8/3i - 1 ) ’ 8 (A) - f + 6/30 ’ 4/3 0 


(3.5) 


where w\ = ^ as shown in (5.f). 


To bound of (3.3), we see it’s important to carefully calculate the numerical flux u „ on 


Proof. 

dK. From the numerical flux formula (3.2) of the DDG schemes, we have hx taken as the element 


diameter or the length of the edge dK. To simplify the proof, here we modify hx to incorporate with 
the mesh geometrical information. We should comment that numerically we observe no difference 
with either choice of hx- 

Again we use edge AB in Figure [l] to illustrate the definition of hx chosen in the numerical 
flux formula (3.2). Let’s have the parametric equation r(t) = fn + xa, 3 , i 6 R to represent the 

And we have points 


normal line through the edge center, 
the normal line with the other two edges of K and A 3 


as the intersection of 
Restricted on edge dK = AB, we take 


and 


hx = h AB = min {||x ^ )3 


— x 


K II) 


|xa' 3 ,3 - x a 3 II} and we have, 


U, 


ds = 


'AB h 


u 


ds + 


Un ds + 


PiK 


ds. 


(3.6) 


JAB JAB ' l ab JAB JAB 

To calculate u nn , we pick two more points along the normal line from each side. We have points 
X A', 5 it = ~h AB ) and x;^(f = — \h AB ) taken in element I\, and points xa ' 3 ,5 (t = h AB ) and x ;<- 3i 4 (t = 
\h AB ) taken in element A' 3 , as shown in Figure [lj Furthermore we denote ux, 1, ■ ■ ■ ,ux, 5 as the 
Uk(x, y) quadratic polynomial solution values on points xa,i, • • •, Xft-,5- As discussed previously, the 
hve points ux, 1 , • ■ ■, ux, 5 are enough to calculate f AB u^ds from the side of element K. Similarly, 
the five points ux 3 ,i, ■ ■ ■ ,ux 3 ,5 are enough to calculate f^u^ds from the side of element A3. 
Essentially the quantity f AB ufrds can be explicitly written out in terms of the ten solution values 
spread out in elements K and A3 as follows, 


I AB h 


u] . fdol AB 


-ds = 


6 h / 


{( u K 3 ,1 + Ux 3 ,2 + ^Ul< 3 , 3 ) — ( U K, 1 + UK, 2 + 4 ua ', 3 )} 


L u " ds = 2 h 


AB {(“3ua 3 ,3 — Ux 3 ,5 + 4«A 3 ,4) + (3ua,3 + Ux, 5 — ^Ux,a)} 

AB 


/ Plh AB [u nn }ds = ' 1 AB {{ux 3 ,3 + Ux 3,5 — 2ua 3 ,4) — (ux ,3 + Ux ,5 — 2ua,4)} • 

JAB U 


(3.7) 


Here we have l AB denoting the length of edge AB. Finally the average uff 1 of (|3.3l) can be written 


out as a function of the solution values that spread out in element A and its neighbors Ai, A 2 , A 3 . 
Thus formally we have, 


—ra+l _ —n 


u 


K 


— u'x H- , / [ Unds + f U n ds + 

area (A) [J AB J BC 


Unds 


'CA 


(3.8) 

















The functional H(---) involves 28 arguments with 15 points from K±, K 2 , K 3 and 13 points from 
A. The first 12 points from A are selected from the calculation of J gK u^ds, see Figure [lj The 


13th one is to be selected by the quadrature rule for cell average see Appendix 5.1 

l K 


Our goal here is to prove solution average G [ m,M] given u^{x,y) G [m,M] on all 

elements. Again, we use a monotone argument showing u 7 ^ 1 is a convex combination of the 
selected solution points. To study the conditions to guarantee that ^(IT, IT, IT, IT) is monotonically 
increasing on the total 28 arguments, it is enough to check out the ten points selected inside A and 


A 3 of (3.8). We first check the five points selected on element A 3 . From (3.7) 

dH 


(3.8) we have, 


dH 


du 


k 3 ,1 


du 


= A 


k 3 ,2 


^AB A) 

h AD 6 : 


dH 


du 


K 3 ,3 


= a H ^-5 + 4A >- 


dH 

duK 3 ,. 


= A 


21 , 


h, 


( 1 — 4/3i), 


dH 

9uk 3) , 


= A 


2 h 


m-i). 


With A = 


At 


area(K) > we on ^A nee( i fio > \ — 6/3i and | < /?i < \ to guarantee the coefficients 
of the five solution values in A 3 being non-negative. Before we carry out the discussion on the 
five points in element K, we need an inequality (refer to Appendix |5.1| ) which reflects geometrical 
property of the mesh partition. 


h, 


= ^tan(min( 6 >i, 0 2 , 0 4 , 6 > 5 )) > ^tan( 0 ). 


From (3.7) and (3.8), we have, 


dH 
9ur, 1 
dH 
9ur, 2 
dH 
du K ,s 
dH 
9uk,A 
dH 


du\ 
du K , 1 
dvfc 
9 u K ,2 

du K ,3 
duK 

9uk,A 

du n K 


dltK, 5 9uk,5 


- A 


A) (l 


6 V h 


\P° 

"6 


- A 


- A 


- A 


l, 


21 , 


2 h, 


l, 


h, 


+ 


+ 




> 


> 


du^ 
du K , 1 
du n K 
9uk/2 


-A 


2/3 0 


-A 


3tan(0) 

2/3o 


=A-5 + 4Aj> At 


3 tan(0) 

4 {Po ~ | + 6/3i) 


- A- 


3 tan(0) 


(1 — 4/3i) > 


m - 1 ) > 


du'x 

9uk,a 

du n K 


-A 


- A 


4(1- 4ft) 
tan( 0 ) 
8/3, - 1 


du K , 5 tan(0) 


To guarantee the coefficients of the five solution points in I\ being non-negative, we need a special 
quadrature rule with all positive weights on all the 12 selected points in A. We refer to Appendix 
5.1 for details of this quadrature rule. The 13th point ua ',13 is selected by the quadrature rule 


also with positive weight. With CFL restriction (3.4)- (3.5) and the quadrature weights (5.1), 
we see H(---) is monotonically increasing on uk, 1 , uk, 2 , uk, 3 , uka and uk, 5 - Similar argument 
applies to edge BC and edge CA, involving solution values in elements Ai and A 2 . Easily we see 
functional H{- ■ ■) is monotonically increasing w.r.t. all 28 point values. With the consistency and 
the monotonicity of H{- ■ ■), we obtain, 

m = H(m, ■■■ , rn) < u ” +1 = H(- ■ ■) < H(M, • • • , M) = M. 


provided that ?n < u''] < (x,y),u'] <i (x,y),u] < 2 (x,y),u , } < 3 (x,y) < M. 


□ 
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3.2 Nonlinear diffusion equation 


In this section we extend the study of 3rd order M-P-S DDG scheme with interface correction (2.3) 
- (2.4) to general nonlinear diffusion equations (2.1) on unstructured triangle mesh. Again, our goal 
is to bound the solution average uff 1 E [ m,M] provided u^(x,y) in the range of [ m,M ]. Take 
test function v = 1 in (2.4) and discretize in time with Euler forward, we have the solution average 
evolving in time as, 

At f , . ,-r~- , , „ At 


Trjn+l _ — n 


U 


K 


= u'k + 


cirea(K) J dK 


(A(u)Vu • n ) ds = u\ + 


area(K) J dK 


Vu • 7 ds. 


(3.9) 



Figure 2: Selected points along direction 7 = A 1 (u^x,;j ))n, for representing the numerical flux 


As discussed previously in (2.3), we apply a new way to calculate numerical flux A(u)Vu • n = 
Vu • 7 with 7 = A T (u) n as a vector pointing from K into its neighbor along the edge. This is true 
because the diffusion matrix A(u) is positive definite and we have 7 ■ n = A T (u) n ■ n > 0 . In a 
word, 7 is a vector always pointing into its neighbor element. To bound uf^ 1 , we need to manage 
to calculate Jq K Vu ■ 7 ds on the three edges of element K. Notice that vector 7 is a nonlinear 
function of the solution, thus we apply a quadrature rule to calculate fg K Vu ■ 7 ds. For example, 
we consider 2-point Gaussian quadrature rule along each edge to approximate the line integral. 
Let’s use point X{j to denote the jth Gaussian point on the ith edge. For each Gaussian point x,;j, 
shown in Figure [2j we see six solution values are enough to calculate the numerical flux, 


u 


Tu • 7|x,. J - W'yl Xij - A) 7 T u*y T Plhi J j\U’y*y\ 


1,3 


(3.10) 


where 7 = A T (u n (xij))ni and hij as the shortest distance from point Xjj to the other edges of 
K and Kj along 7. Again, we bound u r ^ rl by showing it is a convex combination of polynomial 
solution values that spread in element K and its three neighbors K\ , K -2 and K 3 . 

Theorem 3 . 2 . (Nonlinear diffusion equation) Consider DDG scheme with interface correction 


(2.3) - (2.4) with P 2 quadratic approximations on a triangidar mesh. Given u r f-(x,y) in the range 
of we have uf ^ 1 € [m,M] provided, 


A) > 2 “ 




A = 


At 


area(K) 


<A(p 0 ,/3i,e,e). 


(3.11) 
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Again (0o,0i) is the coefficient pair in the numerical flux (2.3), and 6 and 9 as the maximum and 

(3.12) 


minimum angles of the partition Th■ Function A depends on flo, 0i, 6 and 9 as, 

A . 3 -v^ . f 1 1 

3 \ 20o + 80i + 3 80i + 1 

where wq is the minimum quadrature weight in the quadrature rule. 


Proof. Similar to Theorem 3.1 it suffices to study the nronotonicity of uf^ with respect to the 


points selected to evaluate the right hand side of (3.9). Specifically it is enough to study the 


Gaussian points that are used to approximate the line integral on the edges. For one Gaussian 
point Xjj, as shown in Figure [ 2 J six sol ution points selected along 7 = A T (u n (x.ij))ni direction are 
enough to calculate the numerical flux (3.10). We denote u 1 , u 2 , u 3 as the three solution values of 


u^(x,y) selected in K and ul, u 2 , ul as the ones selected in LQ. Again we have h,j denoting the 
shortest distance from point Xjj to the other edges of K and Ki along 7 . 


Now at the Gaussian point Xjj, we can explicitly write out the value of u 7 of (3.10) with 

M = u l - ul 


1 ( 4 u 2 — 3 u 1 — u 3 Au 2 — 3 ul — ul 

U ^ = 2 \ - h - + 


i.j h'Cj ) 2lij_j 


1 


(4u 2 — 3 u 1 — u 3 + 4 u 2 — 3 u\ — u^) 


ui + ~ 2tt ° - +1,3 ~ 21,2 = frK + ul- 2 ul - + » s - 2 » 2 )]. 


hi, 3 
2 


Introduce notation a 1 ’ 3 = 


hj,j 

2 


1,3 


At 


aredijCj^K^j with as the *th edge length and Wj as the jth Gaussian 
point weight. Since we use 2-point Gaussian quadrature rule, the weight ojj = 1. From (3.9), we 
have, 


dul 
dal 
du 2 


(A ,3 duAf 1 


2 a 1 ’ 3 


2 h 


(1 - 4/30, 


•1,3 


a 


1,3 


— ~r (/5o + 40i — -), 


h 


1,3 


Sat 1 

du 3 


a 




dulfl 
du 3 2 h. 

dvJfl 


1,3 


-(8A + 1), 


1,3 


a 


1,3 


du 2 


h. 


1,3 


du 1 


du 1 2 hi. 


-(20 o + 80i + 3). 


With /3o > | — 40i, 4 < 0 , < 1 satisfied in the numerical flux (2.3), we have uffl 1 as a monotone 
increasing function on the solution values u\, uf 0 , ul chosen from element K 3 . Similar discussion 
applies to the solution values used in element K\ and K 2 . With 2-point Gaussian quadrature rule 
approximating the line integral, the quantity uf ^ 1 is monotone increasing with respect to the total 
18 solution values from K \, K 2 and K 3 . 


Similar to Theorem 3.1, we need to have a special quadrature rule for the cell average ufl that 
use all the selected points inside K (18 points in total) with positive weights. We use a similar 
method as the linear case to find such a quadrature rule, see Appendix [5T] Let wq be the minimum 
weight for the selected points. 

With geometry information of the mesh partition, 


/* 

l k 

h 


< 


1 


i,j C sin 9 


C > 


3 — \/3 


(3.13) 
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we see condition (3.11) is sufficient to guarantee the coefficients of solution values used in element 


K being none-negative. Therefore, is monotonically increasing w.r.t. all the selected points 
inside K and its neighbors. Finally we conclude with, 

m < Tfff 1 < M, 


provided that m < u] c (x,y),u] Ci (x,y),u 1 } (2 (x,y),u] <3 (x,y) < M. 


□ 


Implementation of the M-P-S limiter 

Given the quadratic DG polynomial solution u 1 f c (x,y) with cell average u\ 6 [m,M], the 
following limiter ensures u 7 ^(x,y) 6 [■ m,M] for any (x,y) 6 K. 


u\ : ( x, y) = 0{u n K (x, y) - u n K ) + u n K , 


6 = min < 1 


M - u\ 




M k - u n K 

with Mx and mx as the maximum and minimum of u r f^(x, y) on element K. 


m K ~ Uk 


Mk = max u^(x', y), 
{x,y)£K 


rriK = min v^(x,y). 

(x,y)eK 


(3.14) 


(3.15) 


Since Uj^(x, y) is a quadratic polynomial, it is easy to calculate the maximum and minimum value 
over K. 


Notice that the limiter (3.14) does not change the cell average. Moreover, || Uf^(x, y)—u'^(x, y)||oo 
0{h 3 ), if the exact solution is smooth. The proof can be found in [24]. Thus, u] ( (x,y) £ [ m,M] 
has uniform third order accuracy for smooth exact solution. 

Algorithm 3.1. maximum-principle-satisfying DDG scheme with interface correction 


1. At time level t n , we apply M-P-S limiter (3.14) - (3.15) to u 7 j < (x,y) and obtain Uf^(x,y). 


2. Apply DDG scheme with interface correction (2.3) - (2-4) to u r f <: {x,y) and evolve in time with 
SSP Runge-Kutta method to march forward the solution to the next time level t n+ ±. 


For the convection part of (1.1), as discussed in [24j, the solution average at next time level 
is a monotone function with respect to certain solution values (Gauss-Lobatoo points) at time 


level t n . Thus the M-P-S limiter (3.14) - (3.15) with Mk and mx as the maximum and minimum 


of the polynomial solution u r f^(x, y) over the whole element K is enough to guarantee the solution 
average staying in the given bound. 


Remark 3.1. For general convection diffusion equation (1.1), we apply same procedure listed in 


Algorithm 3.1 to guarantee the quadratic polynomial solution stay in the given bound and at the 


same time maintain the 3rd order accuracy. For (1.1), we apply DDG with interface correction 


scheme (2.5) for spatial discretization. 


4 Numerical Examples 

In this section, we present a sequence of examples to demonstrate the performance of M-P-S 
limiter. For all examples in this section, we have the coefficient pair taken as (/3o,/3i) = (5,1) in 


the numerical flux formula (3.2). 
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Example 4.1. Accuracy test on linear diffusion equation 


We start with accuracy check of the DDG with interface correction (3.1) with and without 
M-P-S limiter (3.14) applied on the following linear diffusion equation, 

u t — eAu = 0, (x, y) £ [0,1] X [0,1], t e (0,T), 


with initial data u(x, y, 0) = uq(x, y ) = sin(27r(x + y)) and periodic boundary condition. The exact 
solution is given with, 

u(x, y, t ) = exp(—8 n 2 et) sin(27r(x + y)). 


Here, we take e = 1 and final time t = 0.0001. We implement the scheme with P 2 quadratic 
polynomials on unstructured mesh in Figure |3(a)| and on mesh with obtuse triangles with largest 
angle about r in Figure 3(b) Third order of accuracy is maintained with and without M-P-S 
limiter (3.14) applied, see Table [I] and Table [2| At each time step t n , we set the bounds to be 
= — exp(—8 ir 2 et n ) and ue max = exp(—87 iet n ), which are the minimum and maximum of the 


ue 


exact solution. We use tt m i n and tt max to denote the DG solution minimum and maximum values. 
The overshoots and undershoots are eliminated after the M-P-S limiter applied, see Table [l] and 
Table [2j 



(a) Triangular mesli with h = 0.117. 


(b) Mesh with obtuse triangles, h = 0.148. 


Figure 3: Illustration of meshes 


Example 4.2. Porous medium equation 

In this example we consider the nonlinear porous medium equation 

ut= (« 2 ) xx +(« 2 ) ro , (x,y) G [-10,10] x [-10,10], 


with initial condition 


u(x,y, 0) 


1, (x — 2) 2 + (y + 2) 2 < 6, 
1 , [x + 2) 2 + (y — 2) 2 < 6, 
0, otherwise , 
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h 

without M-P-S limiter 

with M-P-S limiter 


L z error 

Order 

L°° error 

Order 

44min 4ie m in 

Umax 41 e max 

L z error 

Order 

L°° error 

Order 

41min 

l^emin 

Umax 4ie m ax 

0.117 

1.21e-03 


1.10e-02 


-3.91e-03 

3.16e-03 

1.28e-03 


1.10e-02 



0 

0 

0.0587 

1.94e-04 

2.64 

1.28e-03 

3.10 

-2.40e-04 

2.49e-04 

1.95e-04 

2.71 

1.28e-03 

3.10 


0 

0 

0.0293 

2.60e-05 

2.90 

1.59e-04 

3.01 

-1.49e-05 

1.50e-05 

2.60e-05 

2.91 

1.59e-04 

3.01 


0 

0 

0.0147 

3.27e-06 

2.99 

2.01e-05 

2.98 

-9.61e-07 

9.41e-07 

3.27e-06 

2.99 

2.01e-05 

2.98 


0 

0 

0.00733 

4.11e-07 

2.99 

2.52e-06 

2.99 

-5.82e-08 

5.96e-08 

4.11e-07 

2.99 

2.52e-06 

2.99 


0 

0 


Table 1: Accuracy table on triangular mesh of Figure 3(a) 


h 

without M-P-S limiter 

with M-P-S limiter 


L 1 error 

Order 

L°° error 

Order 

4^min 4ie m in 

Umax 41 e max 

Lj error 

Order 

L°° error 

Order 

44min 

^^min 

Umax U&max 

0.148 

8.23e-04 


1.46e-02 


-1.17e-02 

9.70e-03 

1.08e-03 


1.46e-02 



0 

0 

0.0741 

1.23e-04 

2.74 

1.82e-03 

3.00 

-7.86e-04 

5.49e-05 

1.26e-04 

3.11 

2.02e-03 

2.85 


0 

0 

0.0371 

1.68e-05 

2.87 

2.17e-04 

3.06 

-2.66e-05 

4.24e-05 

1.68e-05 

2.90 

2.17e-04 

3.22 


0 

0 

0.0185 

2.14e-06 

2.97 

2.69e-05 

3.02 

-1.12e-06 

9.22e-07 

2.14e-06 

2.97 

2.69e-05 

3.02 


0 

0 

0.00927 

2.70e-07 

2.99 

3.33e-06 

3.01 

-7.17e-08 

1.33e-08 

2.70e-07 

2.99 

3.33e-06 

3.01 


0 

0 


Table 2 : Accuracy table on 


unstructured mesh with obtuse triangles of Figure [3(b)[ 


and zero boundary condition. Piecewise quadratic polynomial solutions implemented on unstruc¬ 
tured mesh (Figure [3 (a) [ ) with size h = 0.00733 are shown in Figure |4j Notice that the minimum of 
the solution is zero. Numerical approximation without M-P-S limiter may become negative which 
may lead the problem ill-posed and cause the computations blow up. Implementations on a coarser 
mesh with h = 0.0587 are carried out, see Figure [5} With M-P-S limiter applied, DDG interface 
correction solutions are maintained strictly inside the bound [ 0 , 1 ]. 

Example 4.3. Strongly degenerate parabolic problem 


We consider the following strongly degenerate parabolic problem with DDG interface correction 
method (2.5), 

ut + ( u 2 ) x + ( u 2 ) y = e{v{u)Vu x ) x + e(v{u)Vu y ) y , (x, y) G [-1.5,1.5] x [-1.5,1.5]. 


Initial condition is given with, 


u(x , y, 0) 


1, (x + 0.5 ) 2 + (y + 0.5 ) 2 <0.16, 
-1, (x-0.5 ) 2 + (y-0.5 ) 2 <0.16, 
0 , otherwise , 


and zero boundary condition is applied. We have e = 0.1 and, 


u(u) 


0, |it| < 0.25, 
1, juj > 0.25. 


Quadratic polynomial implementation with M-P-S limiter is carried out. Here we apply the simple 
slope limiter [3] to compress the oscillations caused from the nonlinear convection term. For the 
slope limiter, we take 7 = 1.5 and M = 5. Implementation on mesh (Figure [3 (a )]) with h = 0.0147 
is shown in Figure [ 6 j The result agrees well those in literature, see mm- 

Example 4.4. Incompressible Navier-Stokes equation in vorticity stream-function formulation 


In this example, we consider to solve two-dimensional incompressible Navier-Stokes equation, 

{ w t + ( uw) x + {vw) y = 

A(f = w, {u,v} = (~<f) y ,(f)x), (4.1) 

(u, v) ■ n = given , (x, y) G 3 D, 
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(c) t = 0.5 


(d) t = 2 


Figure 4: Nonlinear porous medium problem with h = 0.00733. 
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X X 




X 


X 


Figure 5: The cut of the DDG interface correction solution along line x + y = 0atf = 0.005. Red 
circle symbol: no M-P-S limiter. Blue diamond symbol: M-P-S limiter applied. 



Figure 6: Strongly degenerate parabolic problem solution at t = 0.5. 
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written out in the vorticity stream-function format. We focus on the incompressible flow with high 
Reynolds numbers (Re 1), thus explicit treatment on both convection term and diffusion term 
is efficient. 

The initial vorticity is given with w(x,y, 0) = wo(x,y) and periodic boundary condition is 
applied. We have <f>(x,y) denoting the stream function and the velocity field is denoted as (u,v). 
We adopt the method of [TO] by Liu and Shu to solve (4.1). Thus at each time step, we first 


apply P 2 continuous finite element method as the Poisson solver to obtain stream function q i>, then 
have the velocity field and plug them into the vorticity equation and discretize in space with DDG 


interface correction method (2.5), finally update the vorticity DG solution to the next time level. 


High order SSP Runge-Kutta explicit scheme m is applied to march forward solution in time. As 
remarked in [10] that there is a natural match between the vorticity DG solution and the stream 
function. The normal component of velocity field (u, v) ■ n is continuous across all triangle edges, 
thus DG implementation on the convection part of vorticity equation is straight forward. 

We carry out two tests in this example. First one is for accuracy check with exact solution 
maximum and minimum available and being applied with M-P-S limiter at each time step. The 
second one is a vortex patch problem. 

Accuracy Test 


We solve (4.1) with Re = 100. Initial condition is wq(x, y ) = —2 sin(x) sin(y) with D = [0, 27r] x 
[0, 27t] and periodic boundary condition is applied. Exact solution is available with 

w(x,y,t ) = — 2 sin(x) sin(y) exp (— 2t/Re). 


Quadratic P 2 implementations are carried out on mesh Figure 3(a) and on mesh Figure 3(b) with 
obtuse triangles. Errors and orders are listed in Table [3] and Table [4] Again we have we r 


-"max and 

represent exact solution maximum and minimum values. We observe that the M-P-S limiter 
removes all overshoots and undershoots and still maintains the third order accuracy. 


h/27T 

without limiter 

with limiter 


Lr error 

Order 

error 

Order 

1 /; min — Id e min 

W max - We max 

L 1 error 

Order 

error 

Order 

Idmin 

— we m i n 

w max - we max 

0.117 

1.90e-03 


1.85e-02 


-8.17e-03 

7.72e-03 

1.98e-03 


1 .88e-02 



0 

0 

0.0587 

2.68e-04 

2.82 

2.07e-03 

3.16 

-3.64e-04 

6.41e-04 

2.68e-04 

2.88 

2.07e-03 

3.18 


0 

0 

0.0293 

3.62c-05 

2.89 

2.92e-04 

2.82 

-6.40e-05 

-4.50e-06 

3.62e-0-5 

2.89 

2.92e-04 

2.82 


0 

-4.52e-06 

0.0147 

4.54e-06 

3.00 

2.56e-05 

3.51 

-1.29e-06 

3.62e-06 

4.54e-06 

3.00 

2.56e-05 

3.51 


0 

0 

0.00733 

5.84e-07 

2.96 

4.05e-06 

2.66 

-1.53e-07 

1.95e-07 

5.84e-07 

2.96 

4.05e-06 

2.66 


0 

0 


Table 3: Accuracy check on mesh Figure 3(a), final time t = 0.1. 


h/2n 

without limiter 

with limiter 


Lr error 

Order 

error 

Order 

Idmin IPCmin 

IPrnax — "die max 

Lr error 

Order 

L ^ error 

Order 

w 

min Rie m in 

IPmax — tde max 

0.148 

1.31e-03 


2.74e-02 


-8.40e-04 

-7.18e-05 

1.31e-03 


2.74e-02 



0 

-7.18e-05 

0.0741 

1.72e-04 

2.93 

3.35e-03 

3.03 

3.56e-05 

-6.40e-05 

1.72e-04 

2.93 

3.35e-03 

3.03 


4.15e-05 

-6.40e-05 

0.0371 

2.37e-05 

2.86 

4.47e-04 

2.90 

3.36e-06 

1.71e-06 

2.37e-05 

2.86 

4.47e-04 

2.90 


3.36e-06 

0 

0.0185 

2.99e-06 

2.99 

5.33e-0-5 

3.07 

1.79e-07 

1.43e-07 

2.99e-06 

2.99 

5.33e-05 

3.07 


1.79e-07 

0 


Table 4: Accuracy check on 


mesh Figure [3(b)| with obtuse triangles, final time t = 0.1. 


Vortex patch problem 


17 








































Re = 100 

without limiter 

with limiter 

h/2ir 

^Cmin W'e m in 

®max 'fCCmax 

^min ^Cmin ®max 

^^max 

0.117 

-6.96e-01 

6.87e-01 

0 

0 

0.0587 

-2.72e-01 

2 .66e-01 

0 

0 

0.0293 

-8.02e-02 

4.22e-02 

0 

0 

0.0147 

-3.17e-03 

3.00e-03 

0 

0 

0.00733 

-6.41e-04 

7.41e-04 

0 

0 


Table 5: Maximum and minimum of the solutions, Re = 100 at t = 0.1. 


Re = 10000 

without limiter 

with limiter 

Sr- 

to 

=3 

ir m i n We m in 

Wmax. 1Ce ma x 

'ITmin We m i n 

Wmax ^rC m ax 

0.117 

-8.62e-01 

8.62e-01 

0 

0 

0.0587 

-6.14e-01 

7.82e-01 

0 

0 

0.0293 

-5.14e-01 

5.77e-01 

0 

0 

0.0147 

-4.35e-01 

4.79e-01 

0 

0 

0.00733 

-3.52e-01 

2.81e-01 

0 

0 


Table 6: Maximum and minimum of the solutions, Re = 10000 at t = 0.1. 


Now we consider problem (4.1) with initial condition 


wo(x,y) = 


(x,y)€[ §,f]x[|,f], 
(x,y)£[ Ttti 


- 1 , 

1 , 

0 , otherwise , 


L 4 


(4.2) 


and periodic boundary condition. The Reynolds number is chosen to be Re = 100 or Re = 10000. 
We compare the maximum and minimum of the numerical solutions with and without M-P-S 
limiter applied, see Table [5] and Table [6j Mesh Figure 3(a)| is used and quadratic polynomials is 
applied. We also plot the solutions for the case Re = 100 at t = 1, shown in Figure [7j and the 
case Re = 10000 at t = 5, shown in Figure [8j It is clear that the overshoots and undershoots are 
removed with the M-P-S limiter applied. 


5 Appendix 


As discussed in Theorem 3.1, we need to find a 13-point quadrature rule that is exact for P 2 


polynomial and include all the selected points (total 12) as quadrature points. For example we 
consider the edge AB shared by K and K :$. As shown in Figure lj five points x/^* (i = 1, • • • , 5) 
from element K’s side are used in the calculation of (3.6) - (3.7). In this section, our goal is to 


find such a quadrature rule with the selected points as quadrature points and having non-negative 
weights. Specifically we need to find the minimum weight in the quadrature rule, since it is used in 


(3.4) - (3.5) to bound CFL condition and further identify suitable time step size for time evolution. 


For convenience, we introduce notation |A'| to represent the area of triangle element K. 


5.1 Quadrature rules for cell average 

Again we use edge AB to illustrate the way to find such quadrature rule with ~x-K,i (i = 1, ■ • • , 5) 
as quadrature points. The method we investigate is that we look for weights Ui in the following 
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Figure 7: Contours of the solutions, Re = 100, at t = 1 with mesh size h = 0.0293 x 27 t. Top: 
no M-P-S limiter (u^mm = —1.0011, w max = 1.0008); Bottom: add M-P-S limiter (w m i n = —1, 
W maj c = 1). 30 equally spaced contour lines are plotted. 
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Figure 8: Contours of the solutions, Re = 10000 at t = 5 and mesh size h = 0.0293 x 27r. Top: 
no M-P-S limiter (ta m i n = —1.1211, w max = 1.1508); Bottom: add M-P-S limiter (u; m i n = — 1, 
w mSLx = 1). 30 equally spaced contour lines are plotted. 
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Point 

Case 1 

Case 2 

X-K, 4 

1 (2A! 1 \ 

6 [\K\ 3 ) 

> ^ tan(0) cot(0) 


< 2Ai 1 ) 
U*l 3 J 

> ^ tan(0) cot(0) 

X K , 1 

1 • ( Al+A2 Wi 

6 \ A'| Wl 

j tan(0) cot((9) 

1 • ( Ai+a 2 

6 V AT Wl 

j > tan(0) cot(0) 

X K ,2 


j tan(0) cot (0) 

1 

6 ' 

W™ 1 

> "'2 tan(0) cot(0) 

XR, 3 

s • (m Wl ) 

> 

4-3 

0 

0 

+3 

§1“= 

Al 

ir 


> ^ tan(0) cot (8) 

XK, 5 

1 ( 2A 1 +A 2 +A 3 „.. l \ _ w 1 

6 ( | A’| Wl ) ~ 6 

1 ( 2 A 1 +A 0 ... \ — Wl 

6 ( |AT| Wl ) ~ 6 


Table 7: Estimate of the quadrature weights 


format, 



1 

W\ 



5 l 

(■ x , y)dxdy = ^ Wi4(x K ,i) + ^ w*u n K (x* K:j ), 
i= 1 J =1 


with other quadrature points ■ (j = 1 , • • • , i) selected in the following to complete the quadrature 
rule. Notice we will combine these points together with other points from edges BC and CA to 
find the location and weight for the 13-th point. 



Figure 9: Left: Element I\ with its neighbor K 3 ; Right: Two cases of points selected along normal 
vector (Al, A2, A3 denote small triangles area). 


The subsection 5.2 below offers one way to find a quadrature rule on triangle element with 
vortices included as quadrature points (with weight uq in (5.4)). We also use a quadrature rule 
for triangle element only with edge centers with equal weights 1. We divide triangle K into small 
triangles with ~x.K,i (* = 1, ■ • • , 5) being vertices or edge centers and use each rule for a half \v?k- 
In this way, all other vertices and edge centers are selected as quadrature points as well. 

Let’s focus the discussion on point x.k, 4 and xk, 5 > shown in Figure [T] or the blue dots in Figure |9j 
First of all, the location of xra and xr ,5 are determined by elements K and A' 3 . There are two cases: 
(1) min( 0 i, 02 ) > min(# 4 ,(? 5 ), then x ^ 5 is inside K , see Figure [9(b) (2) min(#i, 02 ) < min(# 4 ,# 5 ), 


then x /^5 is on the edge of K, see Figure [9(c) The triangle geometrical information is shown and 
marked in Figure [9j We can estimate the quadrature weights of x^j (i = 1, ■ ■ ■ , 5) in each case, 
see Table 0 
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We apply the same procedure to edge BC and CA to include selected points as quadrature 
points. Collect all data from the three edges, we have the estimate on the weights as follows, 


du\ 


du n K 


dux 


p, , * , o > — tan( 6 »)cot( 6 »), 

OUk, 1 OUK, 2 OUk, 3 6 


f)U n 1 

^->-tan( 0 )cot(e), 


(5.1) 

uuk,o 6 


Now we are ready to find a non-negative quadrature rule with only 13 points inside triangle 
K which include the 12 selected points from numerical flux integral on the edges. Let’s reorder 
the quadrature points, and denote ux.i (* = 1 , • • • , 12 ) as the 12 points selected from evaluating 
f gK u^ds with weights Wi (i = 1, • • • , 12). We also reorder the rest points and denote them as u* K ■ 
( j = 1, • • • , m, where m is an integer), with weights Wj (j = 1, ■ • • , m). Let w 13 = Ej =1 w j> we 

have Ej=i Wj = £}=1 Wj + J2JLi w* = 1. Moreover, ^ Ej=i = Ej=i is a convex 

combination of the points u* K rj (j = 1, • • • , m). By the mean value theory, one can find a point w_jc ,13 

inside the convex hall of the points u* K -,j = 1, • • • , m such that uk.va = Ey =1 ^J^ u *k j • Therefore, 
we have a quadrature rule with positive weights of total 13 points inside triangle K which include 
the 12 selected points: {(ux,j, Wj) \j = 1 , • • • , 12 }. 


5.2 Quadrature rule for triangle element with vertices 

In this section, we design one quadrature rule that is exact for quadratic polynomial P 2 (K) on 
any triangle element K. Especially we include the three vortices and three edge centers in the 
quadrature points set. This work is inspired by the quadrature rule designed in [M] , Specifically 
we are interested in finding the weights before the vertices. 

For convenience, we use the position vectors to denote the three vertices of K: v 1 , v 2 and v 3 . 
Thus, the position vector P of any point P inside triangle K can be described by the barycentric 
coordinates ( 6 , 6 , 6 ), i-e., P = 6 v x + 6 v 2 + 6 v 3 . 

We first consider the quadrature rule on the unit square with vertices coordinates as (—— 3 ), 
(— ^), (^, —|) and (|,);) in the u-v plane, and then we use projections/transformations to map 

it to a quadrature rule on triangle element K. 

Let {u a : a = 1, 2, 3} denote the Gauss-Lobatto quadrature points on [—|, f] with weights w a 
(in Table [ 8 j), which is exact for one variable polynomial of degree 3. We have : j3 = 1,2,3} 
(including the left boundary v 1 = —}) denote the 3-point Gauss-Radau quadrature points on 
[— with weights wp (in Table |9j), which is exact for one variable polynomial of degree 4. For a 
two-variable polynomial p(u,v ), we use tensor product of 3-point Gauss-Lobatto for u and 3-point 
Gauss Radau for v as the quadrature rule on the square. The quadrature points can be written as 


S 2 = { (u a ,vP) : a = 1, 2,3; /3 = 1, 2,3} with weights w a wp, listed in Figure 10(a) 


a 

1 

2 

3 

u a 

1 

2 

0 

1 

2 

— 

1 

— 7T 

1 

^OL 

6 

3 

6 


Table 8 : 3-point Gauss-Lobatto 
quadrature rule 


P 

1 

2 

3 

yd 

1 

2 

to(W6) 

to(1 + V6) 

wp 

1 

9 

^(16 + V 6 ) 

Tj(16->/6) 


Table 9: 3-point Gauss-Radau quadrature rule 


Without loss of generality, we assume the orientation of the three vertices v 1 , v 2 and v 3 is 
marked clockwise. We define the following three functions as projections from the square to triangle 
K, mapping the top edge of the square into one vertex and the other three edges to the edges of 
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K. 


gi{u,v) 

g2 (U,V) 
g3(U,V) 




Under each projection gj (i = 1, 2 or 3), the quadrature points S 2 are mapped onto the triangle K, 
i.e. g i(S 2 ), as in Figure 10(b) - Figure 10(d) . Let S 2 K = gi(5' 2 ) U go(S 2 ) U g3(5 2 ). 




Figure 10: Illustration of quadrature points mapping from rectangle element to triangle element. 


We use g i(i = 1, 2,3) and S 2 to construct our triangle element quadrature rule. Let pk(x, y) be 
a two-variable polynomial of degree 2 with cell average pk = jjU J K pj((x,y)dxdy , then we have, 


1 i 

2 / 2 


Pk = |^[ J K PK( x ,y)dA = j-^| J i J ^PK{.gi{u,v)) 


dgi(u,v) 


d(u,v ) 


dudv, z = 1,2,3 


1 ° ^ ^ ^ 

PK(gi(u, v))2{- - v ) dudv = WaW/3 


2 2 
3 3 3 


a =1 / 3=1 


^^^PK(gi(u a , vf} ))^(l ~ vP)w a W /3 = Y PA-(x)^ x . 


i =1 o=l /3=1 


3 2 


(5.2) 


xSS'lf 


With the three vertices v 1 . v 2 and v 3 orientated clockwise, we have the Jacobian 

2 |iL|(2 — v). Notice that pK{gi{u,v))2{\ — v) is a polynomial of u and v with degree 2 and degree 
3, therefore the quadrature rule on S 2 is exact. 

It is easy to show the weights iv x for quadrature points x 6 are non-negative, then we can 
rewrite (5.2) as a combination of quadrature points, see below, 


PK = Y Pk{k)w x + Ypk(v 1 

x£S 2 k \{vV 2 ,v 3 } i=1 


\Wi. 


(5.3) 


We are interested in the weights {iWj}? l for all three vertices v 1 , v 2 and v 3 . Let us take v 1 for 
example. Notice that g 2 (—— \) and g 3 ( 5 ,— \) are the same point (1,0,0), i.e., v 1 . Therefore, 
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the weight of (1,0,0) is 


w 1 


2 

n 

( 1M 

2 

— 

- - -- 

WlWl + - 

3 

i 2 

V 2 )\ 

3 



W3W1 = - (wi + W 3 ) Wl 


2 

81' 


(5.4) 


Remark 5.1. This section only provide one way to construct a quadrature rule on any triangle, 
with three vertices and edge centers included as quadrature points. The goal is to show that one can 
find such a quadrature rule with positive weights for quadrature points. 
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